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ABSTRACT 


Digital  communication  systems  suffer  from  the  channel  distortion  problem  which  introduces 
errors  due  to  intersymbol  interference.  The  solution  to  this  problem  is  provided  by  equalizers  which 
use  a  training  sequence  to  adapt  to  the  channel.  However  in  many  cases  in  which  a  training  sequence 
is  unfeasible,  the  channel  must  be  adapted  blindly.  Most  of  the  blind  equalization  algorithms  known 
so  far  have  problems  of  convergence  to  local  minima.  Our  intention  is  to  offer  an  alternative 
approach  by  using  extended  Kalman  filtering  and  hidden  Markov  models.  They  seem  to  yield  more 
efficient  algorithms  which  take  the  statistics  of  the  transmitted  sequence  into  consideration.  The 
theoretical  development  of  these  new  algorithms  is  discussed  in  this  thesis.  Also  these  algorithms 
have  been  simulated  under  different  conditions.  The  results  of  simulations  and  comparisons  with 
existing  systems  are  provided.  The  models  for  simulations  are  presented  as  MATLAB  codes. 
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I.  INTRODUCTION 


One  of  the  major  practical  problems  in  digital  communication  systems  is  channel 
distortion  which  causes  errors  due  to  intersymbol  interference.  Since  the  source  signal 
is  in  general  broadband,  the  various  frequency  components  experience  different  steady 
state  amplitude  and  phase  changes  as  they  pass  through  the  channel,  causing  distortion 
in  the  received  message.  This  distortion  translates  into  errors  in  the  received  sequence. 

Our  problem  as  communication  engineers  is  to  restore  the  transmitted  sequence  or, 
equivalently,  to  identify  the  inverse  of  the  channel,  given  the  observed  sequence  at  the 
channel  output.  This  task  is  accomplished  by  adaptive  equalizers. 

Typically,  adaptive  equalizers  used  in  digital  communications  require  an  initial 
training  period,  during  which  a  known  data  sequence  is  transmitted.  A  replica  of  this 
sequence  is  made  available  at  the  receiver  in  proper  synchronism  with  the  transmitter, 
thereby  making  it  possible  for  adjustments  to  be  made  to  the  equalizer  coefficients  in 
accordance  with  the  adaptive  filtering  algorithm  employed  in  the  equalizer  design.  When 
the  training  is  completed,  the  equalizer  is  switched  to  its  decision  directed  mode. 

However,  there  are  practical  situations  where  it  would  be  highly  desirable  for  a 
receiver  to  be  able  to  achieve  complete  adaptation  without  the  cooperation  of  the 
transmitter.  This  can  be  the  case  of  an  unfriendly  receiver.  Also  in  the  communication 
schemes  where  the  channel  parameters  change  with  time,  such  as  in  mobile 
communication  systems,  a  training  sequence  must  be  repeated,  leading  to  waste  in  the 
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channel  utilization.  In  these  cases,  the  equalization  must  be  performed  without  a  training 
sequence.  In  other  words,  the  receiver  is  blind  to  the  specific  transmitted  source 
sequence. 

The  development  of  a  parameter  estimation  algorithm  appropriate  for  a  blind 
equalizer  adaptation  is  hindered  by  the  lack  of  the  reference  signal  by  most  recursive 
schemes.  If  the  source  sequence  values  were  known  to  be  independent  and  identically 
distributed  (i.e.,  i.i.d.)  this  property  could  be  used  to  restore  the  output  of  the  channel 
by  filters  designed  via  estimation  theory  techniques.  The  presumption  is  that  proper 
deconvolution  of  the  received  signal  would  restore  the  source’s  transmitted  sequence. 
To  provide  the  reader  with  a  better  understanding  of  our  approach  to  the  blind 
equalization  problem,  we  will  review  below  the  available  blind  equalization  algorithms 
and  the  problems  with  them. 

A.  AVAILABLE  ALGORITHMS 

1.  Constant  Modulus  Adaptive  Algorithm  (CMA) 

In  the  Constant  Modulus  Algorithm  (CMA),  the  error  between  the  magnitude 
(modulus)  of  the  equalizer  output  and  a  constant  is  recursively  minimized,  resulting  in 
an  algorithm  of  complexity  similar  to  the  Least  Mean  Square  (LMS).  The  motivation 
for  these  methods  is  that  by  restoring  the  modulus  of  the  received  signal,  the  channel 
impulse  response  should  be  implicitly  estimated  and  intersymbol  interference  removed. 
Since  the  magnitude  of  the  equalizer  output  is  independent  of  the  absolute  phase,  the  cost 
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functions  in  these  algorithms  are  independent  of  the  transmitted  data  sequence,  and 
hence,  they  are  capable  of  blind  deconvolving  the  received  sequence. 


The  baseband  model  of  a  digital  communication  system  consists  of  the 
cascade  connection  of  a  linear  communication  channel  and  a  blind  equalizer,  depicted  in 
Figure  (1.1). 


BLIND  !  S(t) 
-4  EQUALIZER  ! _ 

b(t) 

i _ i 

signal 


Figure  1.1  Model  of  Digital  Communication  System 

The  channel  includes  the  combined  effects  of  a  filter  at  the  transmitter,  a 
transmission  medium,  and  a  filter  at  the  receiver.  If  we  assume  it  to  be  linear  and  time 
invariant  or  slowly  varying,  it  is  characterized  by  an  unknown  impulse  response  a  ( t) . 
The  nature  of  the  impulse  response  aft)  is  determined  by  the  type  of  the  modulation 
employed.  We  may  thus  describe  the  input-output  relation  of  the  channel  by  the 
convolution  sum 


S(t) 
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b 
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x(t)=£a,(t)  s(t-i)  V  t 

i 


(1.1) 


where  s(t)  is  the  data  sequence  applied  to  the  channel  input  and  x  (t)  is  the  resulting 
channel  output.  We  assume  that 

Ea,2(fc)  =1  •  (1.2) 

i 

This  implies  the  use  of  an  automatic  gain  control  (AGC)  that  keeps  the  variance  of  the 
channel  output  x  ( t)  constant. 

Let  Jb  (t)  denote  the  impulse  response  of  the  ideal  inverse  filter,  which  is 
related  to  the  impulse  response  a  (t)  of  the  channel  as  follows 

(1.3) 

k 

where  6m  is  the  Kronecker  delta. 


use  an  iterative  deconvolution  procedure  to  complete  an  approximate  inverse  filter 
characterized  by  bi(t) .  The  index  i  refers  to  the  tap-weight  number  in  the  transversal 
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filter  realization  of  the  approximate  filter  [Ref.  1].  Thus  at  the  22th  iteration  we  have  an 
approximately  deconvolved  sequence. 

L 

y( t)  =52  &<( t)  x(t-i)  (1.5) 

im-L 

The  convolution  sum  for  the  ideal  inverse  filter  is  infinite  in  extent  in  that  the  index  i 
ranges  from  -00  to  00.  So  we  may  rewrite  Equation  (1.5)  as  follows; 

y(t)  =^Jb,(t)  x(t-i),  £f(t)-0  for  \  i  \  >L 

1 

or,  equivalently 

y(t)  =£jb,{t)  x ( t-i )  [Jb, ( t)  -£>,]  x(t-i)  .  (1.6) 

i  i 

If  we  let 

v(  t)  =£  [jb,  (t)  -Jb(]  x(t-i),  w;=0  for  \  i  \  >L  (1.7) 

I 

then  we  can  write 

y(t)  =x(  t)  +v(  t)  (1.8) 

where  the  term  v(t)  is  called  "convolutional  noise",  representing  the  residual 
intersymbol  interference  that  results  from  the  use  of  an  approximate  inverse  filter. 

The  inverse  filter  output  y(t)  is  next  applied  to  a  zero  memory  nonlinear 
estimator  producing  estimate  s  ( t)  for  the  datum  s(t)  [Ref.l].  We  may  thus  write 
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§(t)  =  g[y{  t)  ]  . 


(1.9) 


Ordinarily,  we  find  that  the  estimate  s  (t)  is  not  reliable  enough. 
Nevertheless,  we  may  use  it  in  an  adaptive  scheme  to  obtain  a  better  estimate  at  iteration 
t+1.  We  have  a  variety  of  linear  adaptive  filtering  algorithms  to  perform  this 
estimation. 

By  looking  at  the  problem,  we  note  the  following: 

■  The  ith  tap  input  at  the  transversal  filter  at  iteration  t  is  x(t-i) . 

■  Viewing  the  nonlinear  estimate  s  (t)  as  the  desired  response,  and  recognizing 
that  the  corresponding  transversal  filter  output  is  y(t),  we  may  express  the  estimation 
error  for  the  iterative  deconvolution  procedure  as 

e(t)  =  s{t)  -  y(t)  .  (1.10) 

■  The  ith  tap  weight  bi  (t)  at  iteration  t  represents  an  old  parameter  estimate. 

Accordingly,  the  updated  value  of  ith  tap  weight  at  iteration  t+1  is 
computed  as  follows 

Jb,(  t+1)  =Jb,(  t)  +/i  x(t-i)  e(t),  Vi  (1.11) 

where  n  is  the  step  size  parameter.  Equations  (1 .5),(1 .9),(1 . 10)  and  (1.11)  constitute  the 
iterative  deconvolution  algorithm,  known  as  CMA,  for  the  blind  equalization  of  a  real 
baseband  channel.  The  ensemble  averaged  cost  function  corresponding  to  the  tap  weight 
update  Equation  (1.11)  is  defmed  by 
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(1.12) 


Jit)  =E  [  eJ ( t)  ] 

=E  [  (Sit)  -y(t)  )2] 

=E  [  (giyit) )  -yit)  )2]  . 

2.  Decision-Directed  Algorithm 

Figure  (1.2)  presents  a  block  diagram  of  the  equalizer  using  the  decision- 
directed  algorithm.  The  main  difference  between  this  algorithm  and  the  CMA  is  the  type 
of  zero  memory  nonlinearity  imbedded  in  the  equalizer.  Specifically,  the  conditional 
mean  estimation  of  the  blind  equalizer  is  replaced  by  a  threshold  decision  device,  i.e., 

§(  t)  =  dec[y(t)  ]  .  (1.13) 


Figure  1.2  Block  Diagram  of  Decision  Directed  Algorithm 
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For  Binary  Phase  Shift  Keying  (BPSK)  we  may  write 


so 


for  symbol  1 
for  symbol  0 


(1.14) 


dec[y(  t)  ]  =sgn  [y{  t)  ] 


(1.15) 


where  sgn(.)  is  the  signum  function. 

The  Sato  and  Godard  algorithms  [Ref.l],  which  are  generally  used  in 
practice,  are  special  cases  of  the  CMA  and  decision-directed  algorithms.  For  the  Sato 
Algorithm  the  zero  memory  nonlinear  function  can  be  shown  to  be 


s  ( t)  =  y  sgn  [y(  t)  ]  (1.16) 

where 


_  E[s2(t)] 

7  E[|fl(t)  j]  ' 


(1.17) 


For  the  Godard  algorithm  the  zero  memory  nonlinear  function  is  defined  as 


g(fc)=  ly(t')  1  [|y(t)  1*  RP  |y(t)  l'‘,-|y<t)  I*-1]  (i.i8) 

where  p  is  a  positive  integer,  and  Rp  is  a  positive  real  constant  defined  by 


R  =  E[|s(t)  |»] 
P  E[  |  s(  t)  \p] 


(1.19) 
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3.  Higher-Order  Moments 

A  different  approach  is  to  use  higher  order  moments.  Just  to  give  the  general 
idea,  we  know  that  the  data  sequence  s(t)  transmitted  through  the  channel  is  non- 
gaussian.  But  the  output  sequence  x(t) , 

x(  t)  =£  a,  {  t)  s  ( t-i)  V  t 

i 

defined  in  Equation  (1.1),  tends  to  be  Gaussian  if  the  length  of  the  impulse  response  is 
long  enough  (or  at  least  more  Gaussian  than  s  (t)  ).  It  turns  out  that  the  Gaussianity 
of  a  sequence  can  be  measured  by  the  Kurtosis  K  (x  (t) )  associated  with  x  (t)  [Ref.2], 
defined  as 

(1.23) 

K(x(t)  )  =  E  flx(t)  |4  }  -  2  E2  {]x(t)  |2  }  ~  1 25  { x  ( t) 2  » 2  , 

which  uses  the  higher  order  moments  of  x  (t) .  It  can  be  shown  that  K  (x  (t) )  =0  if  and 
only  if  x  (t)  is  Gaussian.  We  can  construct  an  algorithm  by  using  the  Kurtosis  where 
the  impulse  response  of  the  equalizer  b(t)  is  such  that  the  filtered  signal  s  (t)  has  a 
maximum  Kurtosis. 

s(t)  =J2b‘{t)  (1.24) 

i 

The  maximization  is  very  non-linear,  and  more  information  can  be  found  in  [Ref.2]. 
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B.  PROBLEMS  WITH  THE  EXISTING  ALGORITHMS 


1.  Convergence  to  Local  Minima 

One  of  the  major  problems  in  the  algorithms  based  on  CMA,  is  convergence 
to  local  minima.  The  ensemble  averaged  cost  function  is  defmed  by  Equation  (1.12), 

J(t)mE  [  (g(y(t))-y( t))2] 

where  y(t)  denotes  the  received  sequence  defmed  by  Equation  (1 .5).  In  the  linear  case 
of  the  LMS  algorithm,  the  cost  function  is  a  quadratic  (convex)  function  of  the  tap 
weights  and  therefore  has  a  unique  minimum.  By  contrast,  the  cost  function  J(t)  of 
Equation  (1 . 12)  is  a  nonconvex  function  of  the  tap  weights.  This  means  that,  in  general, 
the  error  performance  surface  of  the  iterative  deconvolution  procedure  described  herein 
may  have  more  than  one  local  minimum.  The  nonconvexity  of  the  cost  function  J(t) 
arises  from  the  fact  that  the  estimate  s  (t) ,  is  produced  by  passing  the  linear  combiner 
output  y(t)  through  a  zero  memory  nonlinearity,  and  also  because  of  y(t)  itself  a 
function  of  the  tap  weights. 

Therefore,  nonoptimum  fit  of  the  convergent  equalizer  to  the  channel  inverse 
occurs,  which  lead  us  to  low  performance  of  the  overall  system.  Many  case  studies  have 
been  shown  in  [Ref. 3].  Different  algorithms  have  been  tried,  error  functions  have  been 
derived  and  the  results  have  been  shown  as  plots. 

Also,  if  the  channel  drifts  considerably,  any  decision-directed  algorithm 
looses  track  of  the  channel.  So,  their  usage  is  limited  in  varying  channel  conditions. 
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In  this  thesis  we  address  the  problem  of  applying  optimal  estimation 
techniques  to  the  blind  equalization  problem.  In  Chapter  n  we  introduce  an  algorithm 
that  uses  extended  Kalman  filtering,  in  Chapter  m  we  discuss  an  alternate  approach 
which  requires  a  bank  of  parallel  Kalman  filters,  and  in  Chapter  IV  we  consider  a  new 
algorithm  by  using  the  hidden  Markov  models  as  well  as  the  extended  Kalman  filtering, 
that  combines  the  equalizer  with  the  decoder. 


n.  CMA  LIKE  BUND  EQUALIZATION  ALGORITHMS  BY 
USING  EXTENDED  KALMAN  FILTERING 

As  we  have  discussed  in  the  previous  chapter  blind  equalization  is  an  iterative 
process  mostly  based  on  the  optimal  performance  of  prediction-error-based  recursive 
identifiers.  Since  local  convergence  and  channel  drifting  constitute  a  problem  to  be 
addressed;  the  effectiveness  of  the  algorithms  based  on  CMA,  which  are  known  so  far, 
is  still  object  of  discussion.  These  facts  lead  us  to  consider  a  class  of  blind  equalization 
algorithms  based  on  a  Kalman  filtering  approach.  One  of  the  approaches  is  essentially 
a  variation  of  the  CMA  algorithm  which  has  been  discussed  previously.  Since  the 
Kalman  filter  is  used  for  a  linear  model  of  the  signal,  the  use  of  the  extended  Kalman 
filter  which  takes  nonlinearities  into  account,  and  can  still  be  used  effectively  in  nonlinear 
and  linear  models,  is  considered.  By  this  approach,  we  hope  to  address  the  local 
convergence  problem,  which  is  the  main  concern  for  the  CMA  algorithms.  Before 
presenting  our  approach  we  will  briefly  describe  the  derivation  of  the  Kalman  and 
Extended  Kalman  Filter  equations,  in  order  to  help  to  understand  the  main  idea  behind 
the  new  algorithm. 
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A.  AN  OUTLINE  OF  THE  USE  OF  KALMAN  FILTERS  FOR  STATE 


ESTIMATION  (DISCRETE  TIME  CASE) 

1.  Derivation  of  the  Kalman  Filter  Equations 

We  consider  the  estimation  of  the  states  of  a  system  as  represented  by  a 
system  dynamics  model  of  the  form 

x ( t+1)  =Atx(  t) +Bts  ( t) +w(  t)  (2.1) 

y{t)  =Cex{t) +v{C)  (2.2) 

in  which  w(t)  represents  a  zero  mean  white  noise  disturbance  intensity,  and  v(t)  is 
a  zero  mean  white  noise  corrupting  the  sensor  measurement  signals,  y(t).  Their 
autocorrelations  are  given  by 


E[w(t)  ]  ;E[w(t)  wT(t)  ]  =  Qt  6  (m) 

(2.3) 

E[v(t)]  ;E[v(t)vT(t)]  =  Rt  8  (jn) 

where  Rt  and  Qt  are  positive  definite  matrices  and  6  (m)  is  the  discrete-time  impulse 
function. 


if  m=0 
if  m*  0 


(2.4) 


The  term  s  ( t)  represents  the  system  input,  also  At,  Bt  and  Ct  are  matrices  of 
appropriate  dimension  and  possibly  time  varying. 
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To  develop  a  filter  that  is  recursive  for  real  time  applications  we  define  the 
following  algorithm. 

STEP  1  :  Prediction  of  the  state  estimate  based  on  the  model  and  the 
corresponding  new  measurement,  y(t) . 

St  ( t+1 1 t)  1 1)  +Bts(t)  (2.5) 

STEP  2  :  Update  the  estimation  errors  covariance  matrix,  Pt . 

Pt(t+l|t)  =AtPt(t\t)A?+Qt  (2.6) 

STEP  3  :  Recompute  the  correction  (Kalman)  gain,  Kt. 

JTt(  t  + 1)  =Pe(t  +  l|t)  C*[CtPt(t+l\C)  C*+Rt]  - 1  (2.7) 

Here  we  can  see  that  Kt  is  determined  by  At,  Cv  Qt  and  Rt  only,  not  by  the 
measurement,  y(t) . 

STEP  4  ;  Correction  of  estimate,  based  on  the  error  between  the  new 
measurement  y(t+l)  and  its  prediction,  at  time  t+1. 

^(t+l|ic)  =C^(t+l|t)  (2.8) 

*(t+l|t+l)  =*(t+i|t)  +JCe(t+l)  [y(t+l)  -j?(t+l|t)  ]  (2.9) 

STEP  5  :  Correct  the  estimation  errors  covariance  matrix,  Pt. 

Pt(t+l|t+l)  =  [J-JTe ( t+1)  Ct]Pt( t+l 1 1)  (2.10) 

Also  initial  conditions  must  be  chosen  carefully  where 
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*(0|0)  =E[x(0)] 


(2.11) 


is  the  best  estimate  at  t -0  and, 

P0=cov[x{0)  ]  =E  [bc{0)  -£{0|0)Hx(0)  -^{OlO)}7]  =a\l  (2.12) 

is  a  positive  definite  matrix  generally  in  diagonal  form.  It  shows  the  confidence  on 
initial  conditions. 

2.  The  Extended  Kalman  filter  Equations 

The  Kalman  Filter  was  derived  for  the  case  of  linear  state  and  measurement 
equations.  This  filter  can  be  generalized  to  operate  for  nonlinear  state  equations  and 
measurement  equations  in  a  simple  manner.  This  generalization  will  be  referred  to  as 
the  Extended  Kalman  Filter. 

A  general  non-linear  system  is  in  the  form 

x(t+l)  =f{x{t),s{t),w{t),  t) 
y(t+l)  =f{x{t :)  ,v{t)  ,  t) 

where  all  the  variables  are  consistent  with  the  linear  case.  Then, 

E  [w{  t)  3  ;E  [w(  t)  wT(  t)  ]  =  Qtb  (m) 

E[v(t)]  ;E[v(t)vr(t)]  =  Rt  6  (in) 

are  also  valid  for  this  case.  The  Kalman  filtering  algorithm  can  be  modified  to  account 
for  the  nonlinearities  as  follows: 
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STEP  1  : 


.£{t+l|t)  =f  (.£(t|t)  ,s(t)  ,0,  t) 


SIEEJLi 


pe(t+i|t)  =AtPt{t\t)A?+Qt 


where 


s  _  df{x(t)  ,  s  [  t)  ,w(t)  ,  t) 
*  te(F) 


x(t)  »*(tl  t) 
e!  ri  *s(  t) 
wici  »o 


STEP  3  ; 


JCt  ( t+l)  =Pt  ( t+l  1 t)  ( t+l 1 t)  CZ+BtRcB?)  -1 


where 


c  -  dg{x{t)  ,v(t) ,  t) 
*  dx  ( t) 


|x(t)-*U!t) 
v(t)  -0 


and 


v(  t)  «0 


STEP  4  : 


^(t+l|*)  -sr(*(t+l|t)  ,0,  t) 


*(t+l|t+l)  =*(t+l|t)  +JCt(t+l)  [y( t+l)  -?(t+l|t)  ] 


(2.13) 


(2.14) 


(2.15) 


(2.16) 


(2.17) 


(2.18) 


(2.19) 
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STEP  5  ; 


pt(t+i|t+D  =  [J-irt(t+i)  dt]pt(t+i|t)  v2.20) 

for  the  initial  conditions 

*(0|0)  =E[x(0)  ] 

and 

P0  =  cov[jc(0)  ]  . 

B.  SOLUTION  TO  BLIND  EQUALIZATION  PROBLEM  BY  KALMAN 
FILTERING 

In  this  section  we  show  how  the  extended  Kalman  filtering  can  be  applied  for  blind 
equalization.  Let  xt  be  the  received  signal,  and  Wt  be  the  weights  of  the  equalizer.  The 
purpose  of  the  equalizer  is  to  restore  the  original  information  sequence,  i.e., 

Yt=§t=St 

For  a  binary  phase  shift  keying  (BPSK)  signal  St  €  {  +  1,-1}  and  the  optimal 
equalization  is  achieved  when  /  yt  I  2=1 .  Assuming  the  weights  constant  in  time,  we 
can  write  a  state  space  equation  for  the  observations  as 
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trtU+D  =*t(t) 


(2.21) 


\yt\2=Wt  (t)Xt{C)Xt  (t)  Wt(C)  .  (2.22) 

As  we  can  see,  Equations  (2.21)  and  (2.22)  are  similar  to  Equations  (2. 1)  and  (2.2)  with 
added  nonlinearity. 

By  applying  the  extended  Kalman  fJter  approach  to  the  blind  equalization,  the  new 
algorithm  becomes  as  follows: 

Initial  Conditions: 

*(0|0)  =E[x(0)  ] 

P0  =  COV[x(  0)  ]  =OqJ 


STEP  1  : 


w{t+l\k)  =w{  t) 


(2.23) 


STEP  2  : 


£t(t+l|t)  =lPt{t\t)  IT 


since 


£  _  df{w(t)  ,0,0,  t) 

*  awtl 


=j 


(2.24) 


(2.25) 
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STEP  3  : 


jre(t+i)  =pt(t+ 1 1 1)  £t*,[dtpe(t+i|t)dt*'+d]  (2.26) 

where 

Ct=  dg{W^]{^° =2[wt(t)]  T[xc  { t)  ]  [xc(t)}T  (2.27) 

fi  =  dg(w[t]_L0,t±=Q  (2.28) 

c  8v(t) 

A  small  error  is  needed  for  the  Kalman  filter  in  order  to  work  properly,  so  a  scaler  value 
d  is  added  to  the  Kalman  gain  calculation  at  this  step. 

STEP  4,.: 

92t  ( t+l  1 1)  (t)Xt(  t)  X*  (t)Wt(  t)  (2.29) 

#t(t+l|t+l)  =JV£(t+l|t)  +*t(t+l)  [l-^t(t+l|t)2]  (2.30) 

STEPS  : 

Pt(t+l|t+l)  =  [J-JCe(t+l)dt]Pt(t+l|t)  (2.31) 

C.  FEASIBILITY  OF  THE  BLIND  EQUALIZATION  BY  KALMAN  FILTERING 
As  we  stated  earlier,  an  admissible  source-channel-equalizer  combination  need  not 
always  result  in  an  optimum  fit  of  the  convergent  equalizer  to  the  channel  inverse.  So, 
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we  need  to  study  the  average  behavior  of  the  Kalman  Filtering  Equalizer  algorithm.  In 
this  chapter  we  will  consider  only  linear,  time-invariant,  stable,  channel  models  and 
binary  real  sources,  as  in  [Ref. 3],  The  same  steps  have  been  followed,  in  order  to  make 
a  reasonable  comparison  for  our  new  algorithm. 

Figure  (2.1)  is  an  illustration  of  the  performance  with  various  fixed  equalizer 
parametrizations.  A  zero-mean,  binary  input  sequence  s  (t)  is  applied  to  a  channel. 
The  received  signal  x(t)  is  passed  through  the  equalizer.  To  provide  a  sense  of 
"recent  average"  performance,  the  sequence  of  instantaneous  squared  recovery  errors 
between  s  and  §  is  observed.  An  average  squared  error  of  zero  or  four  denotes  perfect 
performance  for  a  differentially  encoded  signal  and  an  average  squared  error  greater  than 
0.4  but  less  than  3.6  indicates  a  practically  unacceptable  error  rate  of  over  10%. 


Figure  2.1  Smoothed  Squared-Recovery-Error  Time  Histories  for  Different  Fixed  Equalizer 


Parametrizations  [(  b, ,  b, )] 
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In  Figure  (2.1)  all  (b0,b,)  settings  (1  ,-0.8),  (-1,0.8),  (-1,0)  and  (1  ,-0.1)  result  in 
perfect  performance  .  The  squared  recovery  errors  associated  with  (1  ,-0.8)  and  (1  ,-0. 1) 
are  perfectly  zeroed.  The  squared  recovery  error  of  4  associated  with  (-1,0.8)  and  (-1 ,0), 
helps  provide  the  exact  negative  of  s  as  s  at  each  sample  instant.  If  this  figure  is 
compared  with  [Ref.3:Fig.3],  it  can  be  seen  that  the  exponential  rise  does  not  appear  for 
the  Kalman  filtering  algorithm  since  we  did  not  use  a  smoothing  filter.  Also  a  (b0,b,) 
setting  of  (0,1)  can  not  be  seen  in  Figure  (2.1)  since  in  this  case  the  model  is 
undetermined. 

As  a  second  step,  the  performance  of  the  Kalman  filtering  algorithm  for  different 
initializations  is  studied.  The  (b0,b,)  setting  of  (1, -0.6)  is  used  for  the  equalizer  then  the 
initial  values  on  a  circle  of  radius  2  is  applied  to  the  system.  Hie  results  are  shown  in 
Figure  (2.2).  If  compared  to  [Ref.3:Fig.5  and  6]  the  Kalman  filter  does  not  exhibit  the 


Figure  2.2  Parameter  Trajectories  for  Blind  Equalizer  Adaptation  by  Kalman  Filtering  (Initial 
Values  on  the  Circle  of  Radius  2 .Final  Values  in  Black  Boxes,Local  Minima  Are  Small) 
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amount  of  local  minima  shown  by  other  CMA  algorithms.  For  most  of  the  initial 
conditions  the  parameter  vector  converges  to  the  global  minima. 

For  a  perfect  solution  to  the  problem,  we  need  to  understand  the  behavior  of  the 
algorithm  when  convergence  to  local  minima  occurs.  If  we  initialize  (bo.b,)  at  (0,2)  the 
algorithm  converges  to  a  global  minimum.  The  normalized  squared  recovery  errors 
between  s  and  t  for  this  case  show  a  distribution  equivalent  to  N(0,1).  The  situation  is 
shown  in  Figure  (2.3). 
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Figure  2.3  Distribution  of  Squared  Recovery  Errors  When  Algorithm  Globally  Converges 


Also,  we  can  see  from  Figure  (2.2)  that  there  are  initial  conditions  which  lead  to 
local  minima.  As  in  the  example  of  initialization  (0.6, -1.9),  the  algorithm  converges  to 
local  minimum.  In  this  case,  the  normalized  squared  error  between  s  and  $  has  a  larger 
mean  and  variance  than  the  first  case.  Also,  it  does  not  show  a  normal  situation.  The 
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situation  is  shown  in  Figure  (2.4).  This  property  can  be  used  to  determine  the 
convergence  to  a  local  minimum. 


The  studies  up  to  this  step  show  that  blind  equalization  algorithms  using  Kalman 
filtering  seem  to  work  better  than  other  algorithms  in  terms  of  convergence  rate  for 
linear,  time  invariant,  stable  channel  models  and  binary  sources  using  BPSK.  However 
they  still  suffer  from  the  local  convergence  problem  as  seen  for  the  CMA  algorithms. 
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m.  BUND  EQUALIZATION  ALGORITHM  BY  USING 
PARALLEL  KALMAN  FILTERING 

As  seen  in  the  previous  chapter,  the  algorithm  based  on  the  extended  Kalman  filter 
for  blind  equalization,  suffers  from  the  local  convergence  problem,  as  most  of  the  other 
CMA  algorithms.  For  the  algorithm  to  be  applicable,  we  need  to  address  this  problem. 
As  a  candidate  solution  we  choose  a  particular  parallel  Kalman  filtering  approach,  which 
is  described  in  detail  in  [Ref. 4].  In  this  chapter,  we  will  present  this  algorithm  and 
discuss  the  results  of  our  studies  with  this  algorithm. 

A.  DEVELOPMENT  OF  THE  PARALLEL  KALMAN  FILTERING 

ALGORITHM 

This  algorithm  is  an  approximation  to  the  optimum  maximum  a  posteriori  (MAP) 
sequence  estimator  for  a  priori  unknown  channels.  The  sequence  probabilities  can  be 
computed  using  the  innovations  derived  from  a  bank  of  Kalman  filters. 

In  the  discrete-time  channel  and  signal  model,  x(t)  denotes  the  output  of  a 
matched  filter  at  time  t,  a  (t)  is  the  current  transmitted  symbol  and  {an  ( t) }  represent 
the  effective  time  varying  channel  coefficients.  For  BPSK,  the  s  (t)  are  real  valued, 
taking  on  values  {+1  ,-l).  The  channel  coefficients  an  ( t)  represent  the  convolution  of 
the  actual  intersymbol  interference  (ISI)  channel  impulse  response  with  that  of  a 
prewhitening  filter,  which  is  included  to  insure  that  the  additive  noise  samples,  n  ( t) , 


are  uncorrelated. 


(3.1) 


x(t)  =5^a„(t)  sit-n)  +  n(t) 

**0 

The  additive  noise  sequence  is  complex  white  Gaussian  with  variance  oa  2  and  Nb+1  is 
the  length  of  the  channel  impulse  response. 

For  convenience  in  the  derivation,  the  following  sequences  are  also  defmed  for  the 
ith  of  M*+1  possible  sequences,  where  M  is  the  symbol  alphabet  size: 

The  cumulative  measurement  sequence, 

X,={x ( t)  ,x(t-l)  ,  . . .  ,x(0)  } 

A  cumulative  data  sequence, 

D,>,={d1(t),dt.(t- l),...,d,.(0)} 

A  data  subsequence,  comprising  the  data  symbols  associated  with  the  channel 
coefficient  vector  £  (t)  at  time  t, 

d<an,={  d,  ( t)  ,  d,  ( t-1)  ,  . . . ,  d,(  t-Nb)  }  . 

In  the  development  of  the  Kalman  filter  channel  estimator,  it  is  assumed  that  the 
coefficient  bn  (t)  evolve  according  to  the  following  complex  Gaussian 
autoregressive(AR)  process  model 

b(t+l)  =F  b(t)  +w(t)  ,  (3.2) 

with  the  coefficient  vector  defined  by 
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£(t)=[  i>,  ( t)  ,  Jbj  ( t-1) . Jb„ i(0)]r  . 

In  Equation(3.2)  £  is  the  one  step  transition  matrix  and  w(t)  is  a  white  Gaussian 
process  with  covariance  matrix  Q. 

The  optimum  MAP  sequence  estimator  can  be  written  in  the  following  recursive 
form  for  the  assumed  channel  and  signal  models 

p(DjXt)  -1  p(x(t)  | DiJt,X,.x)  (3.3) 

for  i  =  where  p  (Dit  t\Xt)  represents  the  probability  of  the  ith  possible 

data  sequence  given  cumulative  measurements  Xt,  and  c  is  a  normalization  constant. 
The  likelihood  p(x(t)  |  Dit  t,Xc)  is  given  in  terms  of  the  Kalman  filter  innovations 
as 

p(x(t)  |Dfi/,Jrt.,)=N(§I(t)  ,a]{t\ t-1))  ,  (3.4) 

where  N  (x,  ax)  denotes  a  Gaussian  density  with  mean  x  and  covariance  ax2.  The 
estimated  signal  (t) ,  is  given  in  terms  of  the  conditional  channel  estimates  according 
to 

N‘ 

t)  cl  d,(t-n)  ,  (3.5) 

1**0 

where  bi  n(t\t-l)  is  generated  by  the  Kalman  filter  equations  discussed  in  Chapter 

A 

n.  The  estimate  b±  n(t  \  t-l)  can  be  shown  to  be  exactly  equal  to  the  conditional 
mean  of  the  channel  coefficients  an(t),  under  the  AR  process  model  in  Equation(3.2) 
when  conditioned  on  data  sequence  Di  t.  Thus, 
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\  1 1  t~l )  =E  [bH  ( t)  |  DiiltX,.t]  . 


(3.6) 


From  Equations  (3.3)  and  (3.4),  it  is  seen  that  the  optimum  MAP  sequence 
estimator  requires  a  bank  of  M**'  Kalman  filter  channel  estimators,  each  conditioned  on 
a  different  Di  t.  The  MAP  probabilities  of  each  sequence  are  then  obtained  as  a  product 
of  the  corresponding  innovations  likelihoods. 

By  using  the  concept  of  reduced  state  sequence  estimation  (RSEE),  the  algorithm 
becomes  as  follows  (more  detailed  derivation  can  be  found  in  [Ref. 4]): 

i.  Define  observation  vectors 

2i,(t)=[  d,  ( t)  ,  d,  ( t-l ),...,  d,  ( t-Nb)  ]  (3.7) 

ii.  Compute  conditional  innovations  covariances 

<r?(t|  t-l)  =£,( t)  h,T (t)  +  a2n  (3.8) 

iii.  Compute  signal  estimates 

", 

8, it)  =E^(tl t_1)  d.(t*-n)  (3*9) 

n*Q 

iv.  Compute  conditional  measurement  update 

jb^ t|  t)  =bt(t\ t-l)  +— — -  h,T(t)  (x(t)  -§,( t) )  (3.10) 

erf  ( 1 1  t-l) 
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v.  Update  weighing  probabilities 


vi.  Compute  one  step  predictions 


P(Py*».  I*> 


r  p<d.,aix,i 


u.12) 


B.  ADMISSIBILITY  OF  THE  BLIND  EQUALIZATION  BY  PARALLEL 
KALMAN  FILTERING 

In  order  to  be  consistent  with  the  comparison  of  the  algorithms,  we  will  consider 
only  linear,  time  invariant,  stable,  channel  models  and  BPSK.  Then,  s  (t)  e  {+1,-1} 
therefore  M=2  ,and  the  length  of  the  channel  impulse  response  is  Nb+1  =3  (assumed). 
Now  MNb+1=23=8  different  observations  are  needed.  These  are: 

A>(t) 

ij,  (t)  ={-i,  -l,  +1} 
i32(t)={-l,+l,-l} 
th(t)  ={-l,  +1,  +1} 
ii4(t)  ={+i,  -1,-1} 
h j(t)  ={H,-1,  +1} 

23,  (t)  ={+l,+l,-l} 
hj(t)  ={+l,  +1,  +1} 

Assuming  oa2=0.01  (small)  and  a=l  (according  to  [Ref. 4]  0<a<2),  also 
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2no,  ( t|  t-1) 


_  (x(t)  -§,(t)  )\3.i3) 
2ff?(t|  t-1) 


Substitution  of  all  these  parameters  into  the  algorithm  discussed  in  Section  A  yields  the 
performance  with  various  fixed  equalizer  parametrization,  which  is  shown  in  Figure 
(3.1). 


Figure  3.1  Smoothed  Squared-Recovery-Error  Time  Histones  for  Different  Fixed  Equalizer 
Parametrizations  [(b,  ,b,)] 

In  Figure  (3.1)  all  (b^.b,)  settings  result  in  perfect  performance  though  H(z)=£l. 
Also,  the  phase  ambiguity  problem  which  occurs  in  the  extended  Kalman  filtering 
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algorithm  is  solved.  The  exponential  rise  does  not  appear  since  we  did  not  use  a 
smoothing  filter,  and  (b0,bj)  setting  (0,1)  does  not  appear  in  the  figure  since  the  model 
is  undetermined  for  this  case. 


The  performance  of  the  algorithm  for  different  initializations  and  the  (b0,bi)  setting 
of  (1  ,-0.6)  is  also  studied.  The  initial  values  are  shown  as  (x)  and  the  final  values  are 
shown  as  (o)  in  Figure  (3.2).  There  is  no  convergence  to  local  minima  for  the  model, 
so  optimum  results  are  achieved. 


Figure  3.2  Parameters  Initial  and  Final  Values  for  Blind  Equalizer  Adaptation  by  Parallel  Kalman 
Filtering 


However,  during  our  studies  we  observed  that  this  algorithm  does  not  work 
properly  all  the  time.  It  fails  to  converge  to  global  minima  60%  of  the  time.  Also,  it 
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can  not  work  as  well  as  the  extended  Kalman  filtering  approach  in  nonminimum  phase 
systems.  When  the  algorithm  conveiges  to  local  minima,  there  is  no  way  to  determine 
the  situation  with  this  algorithm  as  in  the  Kalman  filtering  approach.  So  the  algorithm 
becomes  unreliable,  since  we  can  not  decide  whether  the  global  minima  have  been 
reached  or  not  by  any  means.  The  usage  of  multiple  Kalman  filters  is  another 
disadvantage  of  this  algorithm,  when  the  cost  is  considered. 
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IV.  APPLICATION  OF  HIDDEN  MARKOV  MODELS  TO  BLIND 
EQUALIZATION  ALGORITHMS  WITH  KALMAN  FILTERING 

In  communication  systems  the  transmitted  sequence  s  (t)  is  not  a  random 
independent  identically  distributed  (i.i.d.)  sequence.  Indeed  forward  error  correction 
(FEC)  coding  is  applied  to  a  message  sequence  u(t) ,  in  order  to  produce  the 
transmitted  sequence  s  (t) .  At  this  point  we  realized  that,  none  of  the  algorithms 
known  so  far  uses  the  statistics  of  the  transmitted  sequence  s  (t) ,  which  might  be  useful 
for  the  channel  equalization.  Obviously,  if  we  decide  to  use  the  statistics  of  the 
transmitted  sequence  s  ( t )  in  equalization,  we  can  combine  the  equalizer  with  the  FEC 
decoder.  It  is  a  fact  that,  most  of  the  FEC  coding  systems  which  are  used  in 
communication  systems,  are  convolutional  coders  and  decoders.  They  mostly  use  the 
Viterbi  soft  decoding  technique,  which  is  an  optimum  estimation  for  Markov  models. 
So  in  our  new  approach  we  will  try  to  combine  the  Kalman  filtering  algorithm  with 
hidden  Markov  models. 

A.  HIDDEN  MARKOV  MODELS 

Consider  a  system  that  may  be  described  at  any  time  as  being  in  one  of  a  set  of  N 
distinct  states  indexed  by  {1,2,3,...,N}.  At  regularly  spaced,  discrete  times,  the  system 
undergoes  a  change  of  state  according  to  a  set  of  probabilities  associated  with  the  state. 
We  denote  the  time  instants  associated  with  state  changes  as  t= 1 ,2,3,... ,  and  we  denote 
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the  actual  state  at  time  t  as  gt.  For  a  discrete  time,  first  order,  Markov  model,  the 
probabilistic  dependence  is  truncated  to  just  the  preceding  state. 


p[q,=j‘l<2r-i=i«  <Z-2=k>  <£- i=1>  • . .]  =p[gf=j|q;-i=i]  <4*1) 

Furthermore,  those  processes  in  which  the  right  hand  side  of  Equation  (4.1)  is 
independent  of  time,  lead  to  the  set  of  state  transition  probabilities  a^j  of  the  form 

a^ptqrjlg.i3-*]  #  Is i,j&N  (4.2) 

with  the  following  properties 


a^O  V  j,i 

*  (4.3) 

E  a«=1  •  v  1 

j‘  i 

The  notation: 


7r(=p[g,=i]  lsisiV  (4.4) 

is  used  to  denote  the  initial  state  probabilities. 

The  concept  of  Markov  models  can  be  extended  to  include  the  case  in  which  the 
observation  is  a  probabilistic  function  of  the  state,  so  the  resulting  model  (which  is  called 
a  hidden  Markov  model)  is  a  doubly  embedded  stochastic  process  with  an  underlying 
stochastic  process  that  is  not  directly  observable  (it  is  hidden)  but  can  be  observed  only 
through  another  set  of  stochastic  processes  that  produce  the  sequence  of  observations. 
As  a  result  a  hidden  Markov  model  (HMM)  for  discrete  symbol  observations  can  be 
characterized  by  the  following  elements: 
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i) .  N,  the  number  of  states  in  the  model, 

ii) .  M,  the  number  of  distinct  observation  symbols  per  state.  We  denote  the 
individual  symbols  as 

v=  {v1,V2,v3, 

iii) .  The  state  transition  probability  matrix  A=  {aij}  where 

Isi,  j*N  (4.5) 

iv) .  The  observation  symbol  probability  distribution,  B=  {bj  (k) },  in  which 

bj(k)  -p[o=vk\q,=j]  IskzM  (4.6) 

defines  the  symbol  distribution  in  state  j,  where  j  =  l,2,3,...,N.  Also  a  typical 
observation  sequence  of  the  model  can  be  shown  as 

°=  (°1  >02,°3,  .  •  .  ,  Ot } 

v) .  The  initial  state  distribution  {■ni}  in  which 

Is  isN  (4.7) 

So  a  complete  specification  of  a  HMM  requires  specification  of  two  model 
parameters,  N  and  M,  specification  of  observation  symbols,  O,  and  the  specification  of 
the  three  sets  of  probability  measures  A,  B  and  it.  The  compact  notation 

X=(A,B,7r)  (4.8) 

is  used  to  indicate  the  complete  parameter  set  of  model.  Given  the  appropriate  values 
of  N,  M,  A,  B  and  it,  the  HMM  can  be  used  as  a  generator  to  give  an  observation 
sequence 
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where  each  observation  oc  is  one  of  the  symbols  from  V,  and  T  is  the  number  of 
observations  in  the  sequence.  This  model  is  the  main  model  used  for  FEC  coding,  and 
it  works  as  follows: 

1.  Choose  an  initial  state  g2=i  according  to  the  initial  state  distribution  ir, 

2.  Set  t=l, 

3.  Choose  ot=vk  according  to  the  symbol  probability  distribution  in  state  i,  i.e. , 
bj(k), 

4.  Transit  to  a  new  state  qt+1=j  according  to  the  state  transition  probability 
distribution  for  state  i,  i.e.,  a±j, 

5.  Set  t=t+l;  return  to  step  3  if  t<T;  otherwise,  terminate  the  procedure. 
Now,  let  us  approach  the  problem  from  our  point  of  view.  We  have  the 

observation  sequence  O,  which  is  produced  by  the  model,  described  above,  and  we  know 
the  model  X.  We  need  to  choose  a  corresponding  state  sequence  q=  (q2 ,  q2,  . .  . ,  qt) 
that  is  optimal  in  some  sense  (i.e.,  best  explains  the  observations). 

We  can  define  the  posteriori  probability  as 


yt(i)  =p[q=i  |  0,  X] 


p[0,  g=i|X] 

pTo|  X) 


p[o,  q=i  |  X] 
£p[0,g=i  I X] 


(4. 10) 


Using  yt(i)  we  can  solve  for  the  individually  most  likely  state  qt*  at  time  t,  as 


q;*  =arg  min, sisJV  [7,(i)  ]  1st &T  (4.11) 

Using  this  criteria  iteratively  the  single  best  state  sequence  (path)  can  be  found. 
To  find  the  single  best  state  sequence  q,  for  the  given  observation  sequence  o,  we  need 
to  define  the  quantity 

(4.12) 

&t  (  i)  *maxfa  #  Q2 9  *  *  • 9  Q- 1 9  Or  9  9 

that  is  the  best  score  (highest  probability)  along  the  single  path,  at  time  t,  which 
accounts  for  the  first  t  observations  and  ends  in  state  i.  By  induction  we  can  compute 
it  recursively  as 

5, ( i )  =  max,  6H(j)  p[q;=i  |  q;.,=i)  p[o,|q;=i]  (4.13) 

To  actually  retrieve  the  state  sequence,  we  need  to  keep  track  of  the  argument  that 
maximized  Equation  (4.13),  for  each  t  and  j.  The  formal  technique  for  doing  the  job 
is  based  on  dynamic  programming  methods,  and  is  called  the  Viterbi  algorithm.  [Ref.5] 


B.  COMBINING  KALMAN  FILTERING  ALGORITHMS  WITH  HMM 

In  general,  a  message  sequence  u  (t)  is  coded  by  a  HMM.  The  block  diagram 
of  an  overall  communication  scheme  is  shown  in  Figure  (4. 1).  The  discrete  coder  states 
z(t)  where  z  (t)  €Z={i,2,3,  — ,L}  and  the  output  sequence  s  (t)  are  dependent 
on  the  older  coder  states  and  the  input  sequence  u  (t)  to  the  coder.  Thus, 


z(t+l)  =  F[z(t)  ,u(t)  ] 
s(t)  =  G[z(t)  ,  u(t)  ] 


(4.14) 
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Figure  4.1  Block  Diagram  of  the  Communication  System 


As  we  stated  in  Chapter  II  before;  the  vector  Xt  represents  the  received  signal  values 
and  the  vector  Wt  represents  the  state  coefficients  of  the  equalizer.  So  we  can  show  the 
output  of  the  equalizer  as 

y(t)  =§(t)  =  W,r(t)  X,(t)+  e(t)  (4.15) 

where  e(t)  is  estimation  error  described  in  Chapter  I. 

In  extended  Kalman  filtering  approach  we  assumed  the  channel  changes  slowly  and 
the  model  is  defined  in  Equation  (2.21).  We  account  for  drifts  in  channel  parameters  by 
the  recursion 
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tff(t+l)  +v(t) 


(4.16) 


Then  combining  the  Equations  (4.14),  (4. IS),  and  (4.16)  we  can  represent  the  overall 
model  of  equalizer  and  decoder  as 

z(t+l)  =  F[z(t)  ,u(t)  ] 
h£(t+l)  -W,(t)  +v(t)  (4.17) 

0=  G[z(t)  ,u(t)]  -  wj{t)  X,(t)+  e(t) 

The  difficulty  with  this  state  space  model  is  that  the  state  [Z(t) ,  W(t)  J  is  a 
mixture  of  a  discrete  component  Z  (t)  and  a  continuous  one  W(t) . 

The  estimation  algorithm  is  based  on  the  fact  that,  given  the  sequences 

Zi(t) = [z(l  ),z(2), . . .  ,z(t)] 

Y,(t)=[y(l),y(2),...,y(t)] 

the  estimate  of  w 

Bl  wt  I  Z,  ,  Yt  ]  (4.18) 

is  well  defined,  and  recursively  computable  by  standard  Kalman  filtering  techniques. 

An  overall  recursion  along  the  same  lines  as  the  estimation  for  hidden  Markov 
models  can  be  devised  for  this  case. 

The  suboptimal  optimization  is  based  on  the  definition 

6,(i)  =  ma>^  p(Z/(t-l)  ,z(t)  =i,  Y,(t) )  (4.19) 

where  ztJ  (t-1)  is  the  optimum  path  up  to  state  z  (t-1)  =j.  By  induction  we  can 
update  6t  (i)  as 
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6,(i)  “max,  Ip(y(  t)  | £/.,  ( t-1) ,  z(  t)  =i,  Y,(  t-1)  )p(z(  t)  =i  |  z(  t-1)  =j)  6,.,  ( j)  ] 

(4.20) 

and  keep  track  of  the  indices  on  the  optimal  path 

j=arg  max  [p(y(  t)  |  Z/_,  ( t-1)  ,  z(  t)  =i,  Yt(  t-1)  )p(z(  t)  =i  |  z(  t-1)  =  j)  6,_,  (j)  ] 

(4.21) 

The  optimal  path  up  to  state  z(t)=i  is  then  updated  as 

%(t)=ZU(t-l)  U  z(t)=i  (4.22) 

Also  we  know  all  the  possible  previous  states  j  and  the  state  coefficient  estimates  of  the 
channel  up  to  this  state  w,.,(j)  is  already  computed  by  Kalman  filtering  as  stated  in 
Equation  (4.18).  So  we  can  determine  the  optimum  estimates  of  the  channel  state 
coefficients  by  the  following  procedure: 

ft,(i|  j)  =^.,(j)  +  [G[j,  u(t)  ]  -  Xt(t)  +  e(t)  (4.23) 

W,(i)=most  likelyj 

C.  FEASIBILITY  OF  HMM  AND  KALMAN  FILTERING  APPROACH  TO 
BUND  EQUALIZATION 

This  new  algorithm  works  almost  perfect  in  linear,  time  invariant,  stable,  channel 
models  and  BPSK.  This  new  algorithm  seems  to  be  more  capable  than  the  previous 
ones,  and  we  applied  it  to  more  complex  channels.  We  have  used  different  feasible  FEC 
encoders  (with  different  rates  and  constraint  lengths),  BPSK  signaling,  additive  white 
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Gaussian  noise  (AWGN)  and  different  signal  to  noise  ratios  (SNR)  in  the  studies.  The 
results  are  presented  in  Figures  (4.2)  through  (4.4). 

Since  the  perfect  equalization  requires 

a(t)  *  b(t)  =  1  , 

the  complete  impulse  response  of  the  combination  of  the  channel  and  the  equalizer  must 
appear  as  an  impulse  function,  which  is  the  case  in  our  studies  (shown  at  the  left  upper 
corner  of  the  figures). 

By  looking  at  the  results,  this  new  algorithm  recovers  the  sent  message  without 
errors  after  a  couple  of  iterations,  which  are  very  small  compared  to  the  length  of  the 
transmitted  sequence.  This  algorithm  works  much  better  in  high  SNR  values.  Errors 
do  occur  if  SNR  drops  under  5  dB.  but  further  studies  are  needed  to  determine  the  exact 
margin.  With  this  kind  of  performance,  this  algorithm  can  also  be  a  candidate  to  address 
the  channel  drifting  problem. 
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Figure  4.2  Channel  and  Message  Estimation  for  Rate=l/2,Con.Len=3  Encoder  with  High  SNR 
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Figure  4.3  Channel  and  Message  Estimation  for  Rate=  1/3,  Con.Len=3  Encoder  with  High  SNR 
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V.  SIMULATION  AND  RESULTS 


We  have  simulated  a  flat  Rayleigh  fading  channel  and  tried  our  new  algorithm  on 
it,  in  order  to  see  its  feasibility  in  overall  communication  systems.  The  block  diagram 
of  our  simulation  scheme  is  shown  in  Figure  (5.1). 


Figure  5.1  Block  Diagram  of  the  Simulation  Scheme  with  DBPSK  and  Flat  Rayleigh  Fading 


In  the  simulation  DBPSK  is  used  as  a  modulation  scheme.  DBPSK  has  hardware 
simplicity  and  does  not  require  phase  coherency.  It  is  popular  for  channels  where  the 
phase  shift  changes  slowly  compared  to  the  bit  duration. 
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Also  the  industrial  standard  FEC  encoder,  with  rate  1/2,  constraint  length  7  and 
(133,171),  connections  is  employed.  The  encoder  and  its  connections  are  illustrated  in 
Figure  (5.2). 


The  flat  Rayleigh  fading  channel  model  is  also  used  and  illustrated  in  Figure  (5.3). 
This  model  corresponds  to  a  single  path  Rayleigh  fading  channel.  In  this  channel  there 
is  no  frequency  selective  distortion,  and  the  channel  is  very  suitable  to  represent  mobile 
communication  schemes.  More  detailed  explanation  about  this  model  and  its 
computational  representation  can  be  found  in  [Ref.  6]. 
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Figure  5.3  Block  Diagram  of  Flat  Rayleigh  Fading  Channel  Model 


The  results  of  the  simulation  using  a  signal  to  noise  ratio  (SNR)  of  13  dB.  is  given 
in  Figure  (5.4).  As  it  can  be  seen  from  the  results  of  our  simulation,  the  combination 
of  equalizer  and  decoder  by  HMM  and  Kalman  filtering  approach  performs  fairly  good 
in  complex  channels  and  under  channel  fading  conditions. 
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VI.  CONCLUSIONS 

An  alternative  algorithm  for  blind  equalization  of  digital  communication  channels, 
is  derived  and  studied  by  using  extended  Kalman  filtering  and  hidden  Markov  models. 
This  new  algorithm  seems  more  efficient  since  it  takes  the  statistics  of  the  transmitted 
sequence  into  consideration.  It  seems  to  be  more  robust  with  respect  to  local  minima, 
and  channel  drifting  problems. 

Hardware  implementation  of  this  algorithm  might  be  cost  effective  since  it 
combines  two  major  components  of  a  receiver,  namely  the  equalizer  and  the  decoder  into 
a  single  component. 

Future  studies  can  determine  the  overall  performance  of  this  algorithm  under 
different  channel  conditions  and  for  different  communication  schemes. 
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APPENDIX  A.  MATLAB  SOURCE  CODES 


The  algorithms  described  in  Chapters  n,m  and  IV  were  implemented  using 
MATLAB  software  and  are  listed  below. 


A,  BUND  EQUALIZER  FOR  LINEAR,  TIME-INVARIANT,  STABLE 
CHANNELS  BY  USING  EXTENDED  KALMAN  FILTER 
Filename  :  eqkalm.m 

Title  :  Blind  equalizer  for  linear, time-invariant, stable  channels  by  using  Kalman  filtering 
algorithm 

Date  of  last  revision  :  15  Jul  1993 

Comments  :  This  program  produces  an  array  of  random  digital  signal  values, 

(either  +  or  -1)  to  represent  the  message  signal  .The  length  of  message  signal 
depends  on  the  sampling  amount  m(of  course  bigger  m  makes  better 
estimation).  Then  this  signal  is  applied  to  a  channel  such  as:[x(t)  x(t-l)  ... 
x(t-n)]  where  the  coefficiants  defined  by  user.To  estimate  the  equalizer 
coefficiants,  a  rectangular  window  over  the  received  signal  values,  is  going  to 
be  used  such  as[x(t)  x(t-l)  ...  x(t-n)]  where  (n+ 1)  represents  the  #  of  values 
taken  into  account  by  the  window  (and  by  the  equalizer). In  this  way  program 
tries  to  predict  real  coefficients  of  channel. If  values  match  equalizer  design 
will  be  perfect. 

Averaged  squared  recovery  error  due  to  the  channel  and  qualizer,performance 
due  to  the  different  initializations  and  normalized  squared  error  distributions 
can  be  plotted. 

Input  variables  : 

f :  The  coefficients  of  the  channel  as  a  vector 
m  :  Sampling  amount 

wt:  Initial  values  of  the  coefficients  at  the  equalizer  as  a  column  vector 
set:  Selection  to  continue  with  the  calculation 
op:  Options  to  plot  different  schemes 
Output  variables  : 

b  :  The  estimated  values  of  the  channel  coefficients  by  equalizer 
Associated  functions  :  None 

clear; 

f=input(’NOW  ENTER  the  coefficients  of  your  channel  as  vector.:’); 

m=input(’  ENTER  the  value  of  m  (sampling  amount) . :’); 

sel=input(’  ENTER  1  to  continue . :’); 
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%%%  Produce  message  %%% 
n«length(f)-l; 
y*=sign(randn(l  ,(m + n))); 

while  sel  *  *  1 

wt«.01*ones((n+l),(m+n)); 

wt(:,n)s=input(’ENTER  the  initial  coefficiants  for  your  equalizer  as  column  vector..:’); 
%%%  Arrangements  %%% 
xt~zeros((n+ 1),1); 
et»zeros(l,m); 
errt=zeros(l,m); 
ay»eye((n+  l),(n+ 1)); 
pt“diag([le-l,le-l]); 
x**zeros(l,(n+l)); 

%  %  %  Pass  message  through  the  channel  %%% 
for  i=n+l:m+n 

x(i)=(-f(2:n+l).*x(i-l:-l:i-n)+y(i))/f(l); 

end 

%%%  Kalman  filtering  %%% 
for  j*n+l:m+n 
xt=(xG:-l:j-n))’; 
wtl*=ay*wt(:  j-1); 
ytl  =wt(:  j-l)’*xt*xt’*wt(:  j-1); 
c=2*wt(:  j-l)’*xt*xt’; 
ptl  =ay*pt*ay’; 
gt=ptl  *c’*inv(c*ptl  *c’ + .01); 
pt=(ay-(gt*c))*ptl; 
wt(:j)=wtl+(gt*(l-ytl)); 
et(:j-n)=l-ytl; 

errt(:  j-n)=et(:  j-n)/sqrt(c*pt  V  +  .01); 
end 

%%%  Results  %%% 
b=wt(:,m+n )’ 
pause 

%%%  Estimation  errors  %%% 
for  k=n+l:n+m 
r(k)*=(b(l:n+l)*x(k:-l:k-n)’); 
s(k)=(y(k)-r(k))"2; 
end 

%%%  Plots  %%% 

disp(’For  the  graphs  you  have  the  following  options:’); 
disp(’ENTER  1  to  see  aver.  squ.  error  due  to  the  chan.  &  equa.’); 
disp(’ENTER  2  to  see  performance  due  to  the  different  initializations’); 
disp(’ENTER  3  to  see  normalized  squared  error  distributions’); 

op=input(’ENTER  0  to  quit . :’); 

while  op  =  =  1 , 
plot(s(n+l:n+m)) 
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grid; 

xlabel(’ Iterations’); 
ylabel(’(u(k)-u(k))  ’); 
pause 

disp(’For  the  graphs  you  have  the  following  options:’); 
disp(’ENTER  2  to  see  performance  due  to  the  different  initializations’); 
disp(’ENTER  3  to  see  normalized  squared  error  distributions’); 

op «input(’ ENTER  0  to  quit . :’); 

end 

while  op*  =2, 
axis([-2  2-2  2]); 
plot(wt(l,:),wt(2,:)) 
grid; 

xlabel(’bO’); 

ylabel(’bl’); 

pause 

disp(’For  the  graphs  you  have  the  following  options:’); 
disp(’ENTER  1  to  see  aver.  squ.  error  due  to  the  chan.  &  equa.’); 
disp(’ENTER  3  to  see  normalized  squared  error  distributions’); 

op=input(’ENTER  0  to  quit . :’); 

end 

while  op==3, 

[n2,x2] =hist(errt,30); 
k2=-4:.01:4; 

gt=  1  /sqrt(2*pi)*exp((-(k2 .  ~2))/2); 
m2=max(gt); 

[xb2,yb2] =bar(x2,((n2/max(n2))*m2)); 

plot(k2,gt,xb2,yb2); 

ort=mean(errt); 

var*(std(errt)A2); 

gtext([’mean  =  ’,num2str(ort)]); 

gtext([’var = ’,num2str(var)]); 

pause 

disp(’For  the  graphs  you  have  the  following  options:’); 
disp(’ENTER  1  to  see  aver.  squ.  error  due  to  the  chan.  &  equa.’); 
dispfENTER  2  to  see  performance  due  to  the  different  initializations’); 


op=input(’ENTER  0  to  quit . :’); 

end 

disp(’  ENTER  0  to  exit  to  matlab . :’); 

sel = input(’  ENTER  1  to  continue . :’); 


end 

end; 


51 


B.  BUND  EQUALIZER  FOR  LINEAR,  TIME-INVARIANT,  STABLE  CHANNELS  BY 
USING  PARALLEL  KALMAN  FILTERS 
%  Filename  :  parkal.m 

%  Title  :  Blind  equalizer  for  linear,time-invariant,stable  channels 
%  by  using  parallel  Kalman  filtering  algorithm 

%  Date  of  last  revision  :  18  Jul  1993 

%  Comments  :  This  program  produces  an  array  of  random  digital 
%  signal  values,  (either  +  or  -1)  to  represent  foe  message 

%  signal. The  length  of  message  signal  depends  on  foe  sampling 

%  amount  m(of  course  bigger  m  makes  better  estimation). Then 

%  this  signal  is  applied  to  a  channel  such  as:[x(t)  x(t-l)  ... 

%  x(t-n)]  where  foe  coefficiants  defined  by  user.To  estimate 

%  foe  equalizer  coefficiants,  a  rectangular  window  over  foe 

%  received  signal  values,  is  going  to  be  used  such  as[x(t) 

%  x(t-l) ...  x(t-n)]  where  (n+ 1)  represents  foe  #  of  values 

%  taken  into  account  by  foe  window  (and  by  foe  equalizer). 

%  Also  each  symbol  in  foe  message  alphabet  is  assigned  to 

%  a  specific  Kalman  filter  .In  this  way  program  tries  to 

%  predict  real  coefficients  of  foe  channel. At  foe  end  correct 

%  channel  parameter  estimations  appear  at  least  at  one  of  foe 

%  Kalman  filters  in  foe  bank.if  values  match  equalizer  design 

%  will  be  perfect.(  Program  works  for  alphabet  size  4  only) 

%  Averaged  squared  recovery  error  due  to  foe  channel  and 

%  equalizer  .performance  due  to  foe  different  initializations 

%  and  normalized  squared  error  distributions  can  be  plotted. 

%  Input  variables  : 

%  f  :  The  coefficients  of  foe  channel  as  a  vector 

%  m  :  Sampling  amount 

%  bi:  Initial  values  of  foe  coefficients  at  foe  equalizer  as 

%  a  column  vector 

%  sel:  Selection  to  continue  with  foe  calculation 

%  sec:  Options  to  plot  different  schemes  for  different  filters 

%  Output  variables  : 

%  b*  :  The  estimated  values  of  foe  channel  coefficients  by 

%  associated  Kalman  filter  *. 

%  Associated  functions  :  None 

clc; 

clear; 

f*=input(’NOW  ENTER  foe  coefficients  for  foe  channel  as  vector(2  elements).:’); 

m=input(’  ENTER  foe  value  of  m  (sampling  amount) . :’); 

sel=*input(’  ENTER  1  to  continue . :’); 

hold  off 
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while  sel=*l 

bi—inputO ENTER  the  initial  coefficiants  for  your  equalizer  as  column  vector..:’); 
n = length  (f)-l; 
y*sign(randn(l  ,(m + n))); 
r»zeros(l,(n+l)); 
zO*ones(l,(n+l)); 
zl»ones(l,(n+l)); 
z2=ones(l,(n+l)); 
z3»ones(l,(n+l)); 
et-zeros(l,m); 
errt=zeros(l,m); 
for  i«n+l:m+n 
r(i)»(f(l:n+l)*y(i:-l:i-n)’); 

%  r(i)*(-f(2:n+l).*r(i-l:-l:i-n)’+y(i))/f(l);  %  Nonrecursive  case 

end 


bO=bi; 
hO=[-l  -1); 
sig0=h0*h0’  +  .01; 
p0=0.25; 

bl=bi; 
hl=[-l  11; 
sigl=hl*hl’  +  .01; 
pi  =0.25; 

b2=bi; 
h2=[l  -1]; 
Sig2=h2*h2’+.01; 
p2=0.25; 


% 

% 

96  FOR  EACH  FILTER 
96 


b3=bi; 
h3=[l  1]; 
sig3=h3*h3’+.01; 
p3=0.25; 


for  j=n+l:m+n 
pOO=pO; 
pi  1  —pi; 

p22=p2; 

p33=p3; 


96  FOR  EACH  FILTER 


sO=sum(bO.*hO’);  % 

bO=bO+((l/sigO)*hO’*(r(j)-sO));  % 

NO=l/(sqrt(2*pi*sigO))*exp((-(rO)-sO).A2)/2*sigO);  96  FOR  EACH  FILTER 
p0=N0*(p00+pll);  96 
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sl»sum(bl.*hr); 

bl-bl+((l/sigl)*hr*(r(j)-sl)); 

N1 « l/(sqrt(2*pi*sigl))*exp((-(r(j)-sl)A2)/2*sigl); 
pl-Nl*(p22+p33); 

s2~sum(b2.*h2’); 

b2  «  b2 + (( 1  /sig2)*h2’  *(r(j)-s2)); 

N2«l/(sqrt(2*pi*sig2))*exp((-(r(j)-s2).A2)/2*sig2); 

p2-N2*(p00+pll); 

s3=sum(b3.*h3’); 

b3  «b3 + ((l/sig3)*h3’*(r(j)-s3)); 

N3  ■  1  /(sqrt(2*pi  *sig3))*exp((-(r0)-s3) .  A2)/2  *sig3); 
p3«N3*(p22+p33); 

pn*pO+pl  +p2+p3;  %  NORMALIZER 

pO-pO/pn;  %FOR  EACH  FILTER 

pl«pl/pn; 

p2«p2/pn; 

p3=p3/pn; 

bO«bO*(pO/(pO+pl))+bl*(pl/(pO+pl));56FOR  EACH  FILTER 
bl  =b2*(p2/(p2+p3))+b3*(p3/(p2+p3)); 
b2=b0; 
b3=bl; 
end 


bO’  %  FOR  EACH  FILTER 

bl’ 

b2’ 

b3’ 

for  k=n+l:m+n 

%  zOCk)=(bO(l:n+l),*r(k:-l:k-n)’);  %  NONRECURSIVE  CASE 

%  zl(k)=(bl(l:n+l)’*r(k:-l:k-n)’);  % 

zO(k) = (zO(k- 1 1  :k-n)*(-b0(2 :  n  + 1 )) + r  (k))/bO(  1 ); 
zl(k)=(zl(k-l:-l:k-n)*(-bl(2.n  +  l))+r(k))/bl(I); 
sO(k)=(l-(zO(k)A2));  %A2; 

seO(k) = sO(k)/sqrt(sigO) ; 
sl(k)=(l-(zl(kr2));  %A2; 

sel  (k)  -  s  1  (k)/sqrt(sigl); 
end 


sec— input(’  ENTER  1  for  1&3;2  for  2&4  channels . :’) 
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while  sec-*l 

*  plot(bO(l,:),bO<2,:),*wo’) 
%  plot(bi(l,:),bi(2,:),’wx’) 

*  axis([-2  2  -2  2]); 

%  hold  on 

%  grid; 

%  xlabel(W); 

%  ylabel(’bl’); 

%  pause 

%  plot(sl(n+l:n+m),’w’) 
%  axis([0  500  -1  5]); 

*  grid; 

%  xlabel(’Iterations’); 

%  ylabel(’(u(k)-u(k))  *); 

%  pause 

[n2,x2] =hist(se0,30); 
k2=-4:.01:4; 


r,% 

% 

% 

%  TRAJECTORIES 
% 

% 

% 

% 

% 

% 

%  RECOVERY  ERRORS 
% 

% 

% 

% 

% 


gt=  l/sqrt(2*pi)*exp((-(k2.A2))/2);  % 
m2=max(gt);  % 

[xb2,yb2]  *=bar(x2,((n2/max(n2))*m2));  *  HISTOGRAMS 
plot(k2,gt,,w-’,xb2,yb2,’w’);  % 
ort=mean(se0);  % 

var=(std(seO)A2);  % 

gtext([’mean=’,num2str(ort)]);  % 
gtext([’var=’,num2str(var)]);  % 
pause  % 

sec=input(’  ENTER  1  for  1&3;2  for  2&4  channels . 

end 


while  sec= —2  %  different  parts  are  as  shown  above 

%  plot(bl(l,:),bl(2,:),’wo’); 

%  plot(bi(l,:),bi(2,:),’wx’) 

%  axis([-2  2  -2  2]); 

%  hold  on 
%  grid; 

%  xlabel(W); 

%  ylabel(’bl’); 

%  pause 

%  plot(sl(n+l:n+m),’w’) 

%  axis([0  500  -1  5]); 

%  grid; 

%  xlabel(Tterations’); 

%  ylabel(’(u(k)-u(k))  ’); 

%  pause 

[n2,x2] =hist(sel  ,30); 
k2=-4:.01:4; 

gt=  l/sqrt(2*pi)*exp((-(k2  A2))/2); 
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^A  Sfi  sfi  \A  vft  vft  vft  vA  \A  ^A  vfi  sA  ^  ^A  \A  sfi  \fi  \A  ^A  ^A  ^A 

CPI  o’!  (y(  ^  ^  ^  ^  “  “  “  "  'j'  O’  ^  O’  ^  ^ 


m2»max(gt); 

[xb2,yb2]  -bar(x2,((n2/max(n2))*m2)); 

plot(it2,gt, ’w-\xb2,yb2, ’w’); 

ort— mean(sel); 

var-(std(sel)A2); 

gtext([’mean> ■’ \num2str(ort)]) ; 

gtcxt([’var  *  \num2str(var)]);  % 

pause 

scc*input(’  ENTER  1  for  1&3;2  for  2&4  channels . :’); 

end 

hold  on 

sel»input(’  ENTER  1  to  continue . :’); 

end 

end; 


C.  BUND  EQUALIZER  FOR  LINEAR,  TIME-INVARIANT,  STABLE  CHANNELS 
BY  USING  KALMAN  FILTER  AND  HMM  ALGORITHM 
Filename  :  try  new.  m 

Title  :  Blind  equalizer  for  linear,time-invariant, stable  channels 
by  using  Kalman  filtering  and  HMM  algorithm 
Date  of  last  revision  :  25  Aug  1993 

Comments  :  This  program  produces  an  array  of  random  digital  NRZ,  BPSK 
signal  values,  (either  +  or  -1)  to  represent  the  message 
signal.The  length  of  message  signal  depends  on  the  sampling 
amount  m(of  course  bigger  m  makes  better  estimation).  Also 
applies  a  FEC  code  by  HMM  rate  choosen  by  the  user  .Then 
this  signal  is  applied  to  a  channel  such  as:[x(t)  x(t-l) ... 
x(t-n)]  where  die  coefficients  defined  by  user  .To  estimate 
the  equalizer  coefficients,  by  Kalman  filtering  and  HMM 
approach. Overall  impulse  response  of  the  system,  sent  and 
recived  sequences  are  plotted  as  well  as  the  errors  due  to 
the  equalizer 
Input  variables  : 

m  :  Length  of  the  message 
Output  variables  :None 
Associated  functions  :  bpskgen.m 
conenc.m 
awgn.m 
eqcondec.m 

m=input(’ ENTER  THE  LENGTH  OF  MESSAGE  YOU  WANT:’); 
mes=bpskgen(m); 
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sig~conenc(roes); 
xn  «reros(size(sig)); 

num»[l, 0.6,0. 4,0.3];  den  =zeros(size(num));  den(l)=l; 
*  num«[l,0];  den«[l,-0.6]; 
xt«filter(num,den,sig);  %  channel 

%xn«xt;  %  no  noise 

xn«awgn(xt, 0.277);  %  13  dB  SNR 

*xn«awgn(xt,0.447);  %  5  dB  SNR 

%xn*awgn(xt,0.577);  %  3  dB  SNR 

%xn»awgn(xt,l);  %  OdBSNR 

h =dimpulse(num,den,20); 

[info,tmax]=eqcondec(xn,h); 


pause 

subplot(223),pIot(mes(l:tmax),,w’),  title(’actual  message’) 
subplot(224),plot(mes(l :tmax)-info(l :tmax),’w’),  title(’errors’) 
end  %  of  program 


function  Y  *  awgn(X, sigma) 

%  AWGN  GENERATOR 

%  07-31-1993 

%  Awgn  is  an  m-file  that  adds  awgn  noise  to  the  matrix  X 
%  Where  sigma  is  standart  deviation  of  the  noise 
%  Ec/No  =  l/(2*sigma*2). 

[rr,cc]=size(X); 

W =randn(rr,cc)+ i*randn(rr,cc); 

Y=X+ sigma.  *W; 

disp(’AWGN  IS  ADDED  TO  SIGNAL’); 


function  u=bpskgen(k) 

%  BPSK  GENERATOR 

%  08-01-1993 

%  This  m-file  accepts  k  the  number  of  bits  that  will  be  returned 
%  in  the  vector  u  which  is  a  BPSK  sequence  of  { + 1  ,-l } 
u=sign(randn(size(l  :k))); 

disp(’A  random  message  is  generated  and  coded  in  bipolar  NRZ  form’); 


function  s=conenc(u) 

%  CONVOLUTIONAL  ENCODER 

%  07-31-1993 

%  This  m-file  is  a  feedforward  convolutional  encoder  for  state  transition 
%  matrix  F,output  sequence  matrix  G  and  input  message  matrix  u.  The  number 
%  of  convolution  schemes  can  be  increased  by  adding  F  &  G  matrices  for  new 
%  rates. 
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disp(’To  use  rates  1/2  convolutional  encoder  ENTER  1  ’); 
disp(To  use  rates  1/3  convolutional  encoder  ENTER  2*); 
disp(To  quit  encoding  ENTER  O’); 

sel»input(”) 
while  sel  =  =  1 , 

F-[l,2; 

3,4; 

1,2; 

3,4]; 

G— [0,-1, 1,-1,  1,1,1,  1,-1; 

0,-l,l,  1,-1,1,-1,-1,  1]; 

x-1; 

s =zeros(size(G(: ,  1)’)); 
for  t»  1  :length(u) 

x  sF(round(x),round(l  .5 + u(t)/2»; 
if  xs  s  1, 

sf sG(: ,  (round  (x) + round(l  .5 + u(t)/2))); 
elseif  xs=2, 

sfsG(:,(round(x)+ 1  +round(l  .5+u(t)/2))); 
elseif  x=  *3, 

sf  =  G(:  ,(round(x) + 2  +  round(l  .5 +u(t)/2))); 
elseif  x==4, 

sf  =  G(:  ,(round(x) + 3 + round(  1 .5 + u(t)/2)»; 
end  %  for  if  statement 
s=[s,sf]; 

end  %  for  for  loop 
dispfMessage  is  encoded’); 
dispfTo  quit  encoding  ENTER  O’); 
selsinput(”) 
end  %  for  while  loop 
while  sel  =  =2, 

F=[l,5; 

1,5; 

2,6; 

2,6; 

3,7; 

3,7; 

4,8; 

4,8]; 

G*  [0,-1, 1,-1,  1,-1,  1,  1,-1,  1,-1,  1,-1; 

0,-l,l,-l,-l,  1,-1,  1,-1,  1,  1,-1,  1; 

0,-1, 1,-1,  1,-1,  1,-1,  1,-1,  1,-1,  1]; 

x-1; 

s  szeros(size(G(; ,  1 )’)); 
for  t=l:length(u) 
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x  -  F(round(x),round(l  .5 + u(t)/2)); 
if  x««l|2, 

sf  *=  G(: ,  (round(x) + round  ( 1.5+ u(t)/2)»; 
el  seif  x=  =3|4, 

sf =G(:,  (round  (x)+ 1  +round(1.5+u(t)/2))); 
elseif  x==5|6, 

sf*G(:  .(round  (x) + 2 + round  ( 1 . 5 + u(t)/2))) ; 
dseif  x==7|8, 

sf  =  G(:  ,(round(x) + 3 + round(  1.5+ u(t)/2)»; 
end  %  for  if  statement 
s=[s,sf]; 

end  %  for  for  loop 
disp(’Message  is  encoded’); 
disp(’To  quit  encoding  ENTER  O’); 
sel*input(”) 
end  %  for  while  loop 


function  [uh,tmax]=eqcondec(x,h) 

%  EQUALIZER  AND  CONVOLUTIONAL  DECODER 

%  07-30-1993 

%  This  m-file  is  a  channel  equalizer  and  convolutional  decoder  for  state 
%  transition  matrix  F, output  sequence  matrix  G  .channel  parameter  matrix 
%  x  and  input  signal  s. 

%  Simply  uses  the  Kalman  filtering  algorithm  with  Hidden  Markov  Models. 

%  The  number  of  convolution  schemes  can  be  increased  by  adding  F  &  G  matrices 
%  for  new  rates 

dispfTo  use  rate  =1/2  convolutional  decoder  ENTER  1’); 
disp(’To  use  rate=l/3  convolutional  decoder  ENTER  2’); 
sel=input(”) 
if  sel  —  =  1 , 

F~[l»2; 

3,4; 

1,2; 

3,4]; 

G= [0,-1, 1,-1,  1,1,-1,  1,-1; 

0,-l,l,  1,-1,1,-1,-1,  1]; 
end  %  for  if  statement 
if  sel  =  =2, 

F*[l,5; 

1,5; 

2,6; 

2,6; 

3,7; 

3,7; 

4,8; 

4,8]; 
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G«[0,-i,i,-i,  1,-1,  l,  l.-i,  1,-1,  i,-i; 

0,-l,l,-l,-l,  l.-i,  l,-l,  l,  l,-l,  l; 

0,-1, 1,-1,  1,-1,  1,-1,  1,-1,  1,-1,  1]; 

end  %  for  if  statement 

nf«length(G(:,l)); 

ns  *  length  (F(:,l)); 

n=  i5;  %  order  of  the  filter 

th=zeros(n,ns); 

th(l  ,:)*B5.0*ones(l  ,ns); 

e2*zeros(ns,l); 

P=l*eye(n); 

tmin = round  (n/nf)  +  1 ; 

tmax=floor(100/nf)-l; 

point=zeros(tmax+  l,ns); 

uhat = zeros  (tmax  + 1  ,ns); 

for  tk=tmin:tmax, 

phi=toeplitz(x(tk*nf+  l:tk*nf+nf),x(tk*nf+  l:-l:tk*nf-n+2))’; 

factl=P*phi; 

fact2=fact1’*phi; 

fact = inv(eye(length(fact2)) + fact2); 

K=factl*fact; 

en2=inf*ones(ns,l);  %  start  error  for  new  layer  at  infinity 


if  ns==4, 
for  j  =  l:ns  %  states 
for  i=  1:2  %  inputs 
m=F(j,i); 

if  m  =  =  1 ,  % 

v=G(:,(m+i))-phi’*th(:j);  % 

elseif  m=  =2,  % 

v=G(:,(m+ 1  +i))-phi’*th(:,j);  %  RATE  1/2 

elseif  m=  =3,  % 

v=G(:,(m+2+i))-phi,*th(:,j);  % 

elseif  m=  =4,  % 

v=G(:,(m+3+i))-phi’*th(:,j);  % 

end  % 


etemp=e2(j)+v,*fact*v; 
if  etemp<en2(m), 
en2(m)=etemp; 
thn(:,m)»th(:j)+K*v; 
point(tk+  l,m)=j; 
uhat(tk+l,m)=i; 
end  %  for  if  statement 
end  %  for  loop  of  i 
end  %  for  loop  of  j 
end  %  for  while  loop 
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if  ns*  *8, 
for  j*l:ns  %  states 
fori*  1:2  %  inputs 
m*F(j,i); 

if  m=  =  112,  % 

v-G(:,(m+i))-phi’*th(:j);  % 

elseif  m==3|4,  % 

v*G(:,(m+  l  +  i))-phi’*th(:  j);  % 

elseif m==5|6,  %  RATE  1/3 

v*G(:,(m+2+i))-phi’*th(:,j);  % 

elseif  m=  =7 1 8,  % 

v=G(:,(m+3  +  i))-phi’*th(:,j);  % 

end  % 


etemp =e2(j) + v’  *fact*v; 
if  etemp  <en2(m), 
en2(m)*  etemp; 
thn(:  ,m)  *th(:  j)+ K*v; 
point(tk+  l,m)=j; 
uhat(tk+  l,m)*i; 
end  %  for  if  statement 
end  %  for  loop  of  i 
end  %  for  loop  of  j 
end  %  for  while  loop 
th=thn; 
e2*en2; 

P  *  P-fact  1  *fact  *fact  1  ’ ; 
end 

[em,m]=min(e2); 

thf=th(:,m); 

%  estimated  message 
while  m  ~  =  0, 
uh(tk)*2*uhat(tk+ 1  ,m)-3; 
m*point(tk+l,m); 
tk*tk-l; 
end 

dispfSignal  is  decoded,  now  look  to  the  graphs’); 
pause 

dg; 

hold  off 
c*conv(h,thf); 

subplot(221),  plot(c,’ow’),title(’impulse  resp  of  ch+eq’) 
%nu*min([length(u),  length(uh)]); 
subplot(222),plot(uh(l:tmax),’w’),  title(’estimated  message’) 
end  %  of  program 
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D.  SIMULATION  OF  THE  HMM  AND  KALMAN  FILTERING  ALGORITHM  IN 


FLAT  RAYLEIGH  FADING  CHANNEL 
%  Filename  :  simul.m 

%  Title  :  Simulation  of  combined  equalizer  and  FEC  decoder  (HMM  and 
%  Kalman  filtering  algorithm)  in  flat  Rayleigh  fading  channel 

%  Date  of  last  revision  :  IS  Sep  1993 

%  Comments  :  This  program  simulates  the  combined  equalizer  and  FEC  decoder 
%  (HMM  and  Kalman  filtering  algorithm)  in  flat  Rayleigh  fading 

%  channel.  Works  similar  to  previuos  programs,  uses  DBPSK 

%  instead  of  BPSK. Industrial  standart,  code  rate  1/2  &  constraint 

%  length  7  FEC  encoder  is  employed 

%  Input  variables  : 

%  m  :  Length  of  the  message  sequence 

%  Output  variables  :None 
%  Associated  functions  :  awgn.m 
%  bpskgen.m 

%  conenc64.m 

%  difecod.m 

%  pskmod.m 

%  fade.m 

%  pskdmod.m 

%  eqcod64.m 

m = input(’ENTER  THE  LENGTH  OF  MESSAGE  YOU  WANT:’); 

mes=bpskgen(m); 

erne = conenc64(mes ,  F64 ,  G64) ; 

len=length(cme); 

dcm=difecod(cme); 

sig =pskmod(dcm,7); 

sig  1=  sig’; 

sig2=sigl(:); 

sig3 =exp(i*pi/4)*sig2; 

disp(’POWER  IS  SPLITTED  TO  I  &  Q  CHANNELS’); 
I=fade(0.01,(len+1)*7); 

Q=fade(0.01,(len+ 1)*7); 
sig4 = I.  *real(sig3)  +  j  *(Q.  *imag(sig3)); 
dispfRayleigh  fading  is  applied  to  signal’); 
sig5 =awgn(sig4 ,0.277); 

[emm,  a] = pskdmod  (sig  ,7,len) ; 

[info,tmax]=eqcod64(cmm,F64,G64); 

pause 

subplot(211),plot(mes(l:tmax),’w’),  title(’actual  message’) 
subplot(212),plot(info(l  :tmax),’w’),title(’ received  message’) 
pause 

print  -dmeta; 
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clg; 

subplot(21  l),plot(mes(l  :tmax)-info(l  :tmax),’w’),  titlef'errors’) 
end  %  of  program 


function  s=conenc64(u,F64,G64) 

%  CONVOLUTIONAL  ENCODER 

%  08-30-1993 

%  This  m-file  is  a  feedforward  industrial  stand  art  convolutional  encoder 
%  for  state  transition  matrix  F64,output  sequence  matrix  G64  and  input 
%  message  matrix  u.  F64  &  G64  matrices  are  given  as  a  file  named  incod 64. mat 
%  Be  sure  to  load  that  file  before  working  with  this  function 

X  es  1  * 

s=zeros(size(G64(:,  1)’)); 
for  t=l:length(u) 

x = F64(round(x),round(l  .5 + u(t)/2)); 
if  x=  =  l  1 9, 

sf =  G64(: ,  (0 + round(  1 .5 + u(t)/2))); 
elseif  x==2|  10, 

sf=G64(:,(2+round(l  .5+u(t)/2))); 
elseif  x==3|  11, 

sf*G64(:,(4+round(l  .5+u(t)/2))); 
elseif  x=  =4 1 12, 

sf=G64(:,(6+round(l  .5+u(t)/2))); 
elseif  x===5|  13, 

sf=G64(:,(8+round(l  .5+u(t)/2»); 
elseif  x==6|  14, 

sf=G64(:,(10+round(l  ,5+u(t)/2))); 
elseif  x==7|  15, 

sf=G64(:,(12+round(l  .5+u(t)/2))); 
elseif  x=  =8 1 16, 

sf=G64(:,(14+round(l  .5 +u(t)/2)»; 
elseif  x=  =  17|25, 
sf=G64(:,(16+round(1.5+u(t)/2))); 
elseif  x=  =  18j26, 
sf=G64(:,(18+round(l  .5+u(t)/2))); 
elseif  x=  =  19 1 27, 
sf— G64(:  ,(20 + round(  1 .5 + u(t)/2))); 
elseif  x==20|28, 
sf=G64(:,(22+round(l  .5+u(t)/2))); 
elseif  x=  =21 129, 
sf = G64(:  ,(24 + round(  1 ,5+u(t)/2))); 
elseif  x==22|30, 
sf  =G64(:,  (26+ round  (1 .5+u(t)/2))); 
elseif  x=  =23 1 31, 
sf =G64(:  ,(28 + round(l  .5 + u(t)/2))); 
elseif  x==24|32, 


63 


sf»G64(:,(30+round(l  ,5+u(t)/2))); 
elseif  x==33|41, 
sf=G64(:,(32+round(l  .5+u(t)/2))); 
clseif  x==34|42, 
sf=G64(:,(34+round(l  ,5+u(t)/2))); 
el  seif  x=  =35 1 43, 
sf =G64(:  ,(36+ round(l  .5 +u(t)/2))); 
elseif  x=  =36|44, 
sf  *  G64(:  ,(38 + round(l  .5 + u(t)/2))); 
elseif  x=  =37 145, 
sf = G64(:  ,(40 + round(  1 .5+u(t)/2))); 
elseif  x=  =38 146, 
sf =  G64(:  ,(42 + round(l  .5 + u(t)/2))); 
elseif  x==39|47, 
sf  =  G64(:  ,(44 + round(l  ,5+u(t)/2))); 
elseif  x==40|48, 
sf=G64(:,(46+round(l  .5+u(t)/2))); 
elseif  x=  =49 1 57, 
sf  =  G64(:  ,(48 + round  ( 1 .5+u(t)/2))); 
elseif  x==50|58, 
sf  =  G64(: ,  (50 + round  ( 1 . 5 + u(t)/2))); 
elseif  x=  =51 159, 
sf=G64(:,(52+round(l  ,5+u(t)/2))); 
elseif  x=  =52 160, 
sf=G64(:,(54+round(l  .5+u(t)/2))); 
elseif  x==53|61, 
sf  =  G64(:  ,(56 + round(  1 .5+u(t)/2))); 
elseif  x==54|62, 
sf  =  G64(:  ,(58 + round(l  ,5+u(t)/2))>; 
elseif  x==55|63, 
sf=G64(:,(60+round(l  ,5+u(t)/2))); 
elseif  x==56|64, 
sf  =  G64(:  ,(62 + round(l  .5+u(t)/2))); 
end  %  for  if  statement 
s=[s,sf]; 

end  %  for  for  loop 

disp(’FEC  is  applied  to  Message(ENCODED)’); 


function  [uh,tmax]=eqcod64(x,F64,G64) 

%  EQUALIZER  AND  CONVOLUTIONAL  DECODER 

%  08-2M993 

%  This  m-file  is  a  channel  equalizer  and  convolutional  decoder  for  state 
%  transition  matrix  F64,output  sequence  matrix  G64, channel  parameter  matrix 
%  x  and  input  signal  s. 

%  Industrial  standart  Rate  1/2  code  with  constraint  length  7  is  used  F64  & 

%  G64  matrices  are  defined  in  incod 64. mat.  Be  sure  it  is  loaded. 
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%  Simply  uses  the  Kalman  filtering  algorithm  with  Hidden  Markov  Models. 

%  The  number  of  convolution  schemes  can  be  increased  by  adding  F  &  G  matrices 

%  for  new  rates 

nf  ■  length(G64(: ,  1 )); 

ns«length(F64(:,l)); 

n=  IS;  %  order  of  the  filter 

th*zeros(n,ns); 

th(  1 , :)  «5.0*ones(  1  ,ns); 

e2=*zeros(ns,l); 

P=l*eye(n); 
tmin= round  (n/nf)  + 1 ; 
tmax  *  floor(100/nf)- 1 ; 
point =zeros(tmax+  l,ns); 
uh  at = zeros  (tmax  + 1  ,ns); 
for  tk=tmin:tmax, 

phi*toeplitz(x(tk*nf+  l:tk*nf+nf),x(tk*nf+  l:-l:tk*nf-n+2))’; 

factl  *P*phi; 

fact2=sfactl’*phi; 

fact = inv(eye(length(fact2)) + fact2) ; 

K=factl*fact; 

en2 = inf*ones(ns ,  1 ) ;  %  start  error  for  new  layer  at  infinity 

for  j=l:ns  %  states 
for  i=  1:2  %  inputs 
m=F64(j,i); 
if  m=  =  1 19, 

v=G64(:,(0+i))-phi’*th(:j); 
elseif  m=  =2 1 10, 
v=G64(:,(2+i))-phi’*th(:j); 
elseif  m==3|  11, 
v=G64(:,(4+i))-phi’*th(:  j); 
elseif  m==4|  12, 
v=G64(:,(6+i))-phi’*th(:j); 
elseif  m==5|  13, 
v  «  G64(:  ,(8 + i))-phi  ’  *th(:  ,j); 
elseif  m==6j  14, 
v=G64(:,(10+i))-phi’*th(:,j); 
elseif  m=  =7|  15, 
v*=G64(:,(12+i))-phi’*th(:,j); 
elseif  m==8]  16, 
v=G64(:,(14+i))-phi’*th(:j); 
elseif  m=  =  17|25, 
v=G64(:,(16+i))-phi’*th(:,j); 
elseif  m =  =  1 8 1 26, 
v=G64(:,(18+i))-phi’*th(:J); 
elseif  m=  =  19 1 27, 
v=G64(:,(20+i))-phi’*th(:  j); 
elseif  m==20|28, 
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v«G64(:,(22+i))-phi’*th(:j); 
elseif  m««21}29, 
v=G64(:,(24+i))-phi**th(:j); 
elseif  m==22|30, 
v=G64(:,(26+i))-phi,*th(:  j); 
elseif  m==23|31, 
v=G64(:,(28+i))-phi’*th(:j); 
elseif  m*** 24 1 32, 
v=G64(:,(30+i))-phi’*th(:j); 
elseif  m=  =33 1 41, 
v=G64(:,(32+i))-phi’*th(:j); 
elseif  m= =34)42, 
v = G64(:  ,(34 + i))-phi  ’*th(:  j); 
elseif  m=  =35)43, 
v=G64(:,(36+i))-phi’*th(:  j); 
elseif  m=  =36)44, 
v=G64(:,(38+i))-phi”*th(:  j); 
elseif  m= =37  (45, 
v=G64(:,(40+i))-phi’*th(:j); 
elseif  m==38|46, 
v=G64(:,(42+i))-phi**th(:j); 
elseif  m== 39 1 47, 
v=G64(:  ,(44 + i))-phi’*th(:  ,j); 
elseif  m==40|48, 
v =G64(:  ,(46+ i))-phi’  *th(:  ,j); 
elseif  m==49)57, 
v=G64(:,(48+i))-phi’*th(:,j); 
elseif  m=  =50)58, 
v=G64(:,(50+i))-phi’*th(:j); 
elseif  m==51 159, 
v=G64(:,(52+i))-phi’*th(:,j); 
elseif  m==52|60, 
v = G64(:  ,(54 + i))-phi  ’  *th(:  ,j); 
elseif  m=  =53)61, 
v=G64(:,(56+i))-phi’*th(:,j); 
elseif  m==54|62, 
v=G64(:,(58+i))-phi’*th(:,j); 
elseif  m=  =55)63, 
v=G64(:,(60+i))-phi’*th(:  j); 
elseif  m=  =56|64, 
v  =  G64(:  ,(62 + i))-phi  ’  *th(:  ,j); 
end 

etemp =e2(j)+v’*fact*v; 
if  etemp  <en2(m), 
en2(m)= etemp; 
thn(:  ,m)  =th(:  j) + K*v; 
point(tk+l,m)=j; 
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uhat(tk+l,m)*i; 
end  %  for  if  statement 
end  %  for  loop  of  i 
end  %  for  loop  of  j 
th=thn; 
e2*en2; 

P  *  P-fact  1  *fact*ftct  1’ ; 
end 

[em,m]=min(e2); 

thf«th(:,m); 

%  estimated  message 
while  m  ~  =  0, 
uh(tk)=2*uhat(tk  + 1  ,m)-3; 
m*point(tk+l,m); 
tk=tk-l; 
end 

disp(’Signal  is  decoded’); 


function  out=pskmod(in,N) 

%  PSK  MODULATOR 

%  09-03-1993 

%  This  M  file  accepts  the  differentially  encoded  vector  in 
%  and  #  of  samples  per  symbol  N.  Then  generates  a  matrix 
%  with  dimensions  (length(in)xN) 
out-in’*ones(l,N); 

dispfMESSAGE  DPSK  MODULATED  ); 


function  [out,a]=pskdmod(sig,N,m) 

%  PSK  DEMODULATOR 

%  09-02-1993 

%  This  function  accepts  fading  and  AWGN  effected  sig, 

%  #  of  samples  per  symbol  N  and  total  #  of  symbols  m 
%  Then  for  each  symbol  takes  the  difference  and  sum  of 
%  present  and  past  symbol  sums  and  squares  their  absolute 
%  values.  If  difference < sum  decides  input  bit  is  0(rep.  -1) 
o— ones(l,m); 
for  i=2:m+l, 

dif(i- 1 ) = abs(sum(sig(i ,  :))-sum(sig(i- 1 ,:))).  A2 ; 
su(i-l)=abs(sum(sig(i,:))+sum(sig(i-l ,:))).  A2; 
met(i- 1 ) = dif(i- 1  )-su(i- 1 ); 
if  met(i-l)<0; 

o(i-l) — 1; 
end 
end 

a=met; 
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out»o; 

dispCSIGNAL  IS  DEMODULATED’); 


function  difout=difecod(n) 

*  DIFFERENTIAL  ENCODER 

%  OS-31-1993 

%  This  M  file  accepts  an  input  vector  n(which  is  bipolar  NRZ 
%  coded  form  of  output  sequence)  and  differentially  encodes  it 
%  with  a  reference  bit  1  (which  will  be  inserted  at  the  beginning 
%  of  the  vector, 
y  *  [1  ,zeros(l  ,length(n))]; 
for  i“l:length(n), 
if  n(i)===-l, 
n(i)=0; 
end 

y(i+  l)=xor(n(i),y(i)); 
end 

for  j=l:length(n), 
if  n(i)=  =0, 

nO)— i; 

end 

ify(i)==0, 

yCO— 1; 

end 

end 

difout=y; 

disp(’MESSAGE  DIFFERENTIALLY  ENCODED’); 


function  out— fade(dps,n) 

%  FADING  GENERATOR 

*  09-04-1993 

%  This  function  accepts  differential  phase  shift  dps  and 

%  #  of  samples  n.  Then  creates  n  normal  distributed 

%  variables  passes  these  through  two  RC  filters 

s=exp(-dps/2. 146193); 

ss— ((l-sA2)A3/(l  +sA2))A.25; 

b=[ss]; 

a=[l  -s]; 

t=randn(n,l); 

k=filterO>,a,t); 

out=filter(b,a,k); 
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